Fonte: Storymap UFRRJ

Introdução à Análise Estatística Espacial

O que é Análise Estatística Espacial ?

  • São métodos estatísticos que levam em consideração a localização espacial do fenômeno estudado;

  • Segundo Bailey & Gatrell (1995), “Análise estatística espacial pode ser realizada quando os dados são espacialmente localizados e se considera explicitamente a possível importância de seu arranjo espacial na análise ou interpretação dos resultados”

  • A primeira pergunta a ser feita é: A distribuição dos dados apresenta um padrão aleatório ou apresenta uma agregação definida (clusters) ?

Origem da Estatística Espacial

  • Dr. John Snow (1813-1858) \(\rightarrow\) Considerado pai da Epidemiologia Moderna

Mapeamento dos casos de coléra (\(\bullet\)) e as bombas de água (X) em Londres, 1854.

Objetivos da Estatística Espacial

  1. Análise de padrões espaciais e espaço-temporais (ex: Análise Exploratória Espacial, Correlação Espacial)

  2. Modelagem Espacial (ex: Modelos Estatísticos de Regressão Espacial)

Dependência Espacial ou Autocorrelação Espacial

  • “Independência é um pressuposto muito conveniente que faz grande parte da teoria estatística matemática tratável. Entretanto, modelos que envolvem dependência estatística são frequentemente mais realísticos. . . . Neste tipos de dados a dependência está presente em todas as direções e fica mais fraca a medida em que aumenta a dispersão na localização dos dados.” (Cressie N.,1991)

  • “Todas as coisas se parecem, coisas mais próximas são mais parecidas que aquelas mais distantes.” (Tobler, 1979). Também conhecida como \(1^a\) Lei da Geografia

Tipologia dos Dados Espaciais

Os diferentes tipos de dados espaciais são tradicionalmente classificados de acordo com uma tipologia. Esta caracterização diz respeito a natureza estocástica da observação.

Noel Cressie (1993) divide a estatı́stica espacial em 3 grandes áreas:

  1. Dados de processos pontuais

  2. Dados de geoestatística

  3. Dados de área

Existem métodos estatı́sticos diferentes para descrever ou analisar estes tipos de dados. Exemplo:

Padrões Pontuais

  • O principal interesse está no conjunto de coordenadas geográficas representando as localizações exatas de eventos.

  • Exemplos: Localização de crimes, localização da residência dos casos de dengue, localização de espécies vegetais, etc.

  • Neste caso, o dado aleatório de interesse é a localização espacial do evento.

Estimativas de Kernel (ou Mapas de Calor)

  • Estimador de intensidade de distribuição de pontos:

\[\hat{\lambda}_{\tau}(u) = \dfrac{1}{\tau^2}\sum k(\dfrac{d(u_i , u)}{\tau}), d(u_i , u) \leq \tau\] *Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

  • Localização da ocorrência de casos de dengue em Belo Horizonte/MG

Geoestatística

  • Podemos definir como sendo uma análise de um atributo espacialmente contínuo amostrado em localizações fixas.

  • Os dados compreendem um conjunto de localizações (em geral latitudes e longitudes), mas atribuidos a eles uma medida, como por exemplo: volume de chuva, concentração de poluentes no ar, número de ovos de Aedes postado em ovitrampas, etc.

  • A determinação experimental do semivariograma, para cada valor de \(x_i\), considera todos os pares de amostras \(z(x)\) e \(z(x+h)\), separadas pelo vetor distância \(h\), a partir da equação:

\[\hat{\gamma}(h) = \dfrac{1}{2N(h)}\sum^{N(h)}_{i=1}[z(x_i) - z(x_i + h)]^2 \]

  • Sendo \(\hat{\gamma}(h)\) o semivariograma estimado e \(N(h)\) o número de pares de valores medidos, \(z(x)\) e \(z(x+h)\), separados pelo vetor \(h\).

  • Esta fórmula, entretanto, não é robusta. Podem existir situações em que variabilidade local não é constante e se modifica ao longo da área de estudo (heteroscedasticidade).

*Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

  • Exemplo: Mapa sobre o teor de argila no solo.

Dados de Área

  • Na análise de áreas o atributo estudado é em geral resultado de uma medida sumáraria (ex: contagem, taxas, médias, etc.) em uma determinada área bem definida, por exemplo um polígono.

  • Tal polígono cuja forma pode ser complexa bem como as relações de vizinhança.

  • O objetivo é a detecção e explicação de padrões e tendências observados entre áreas.

Mapa Temático

  • Taxas de câncer de pulmão na população branca masculina nos Estados Unidos, por condados no ano de 1998

Matriz de Vizinhança

Moran Global, Moran Local e Lisa Map

Moran Global

Moran Global

Moran Local

Moran Local

  • Desigualdade no nível distrital na cobertura de saúde reprodutiva, materna, neonatal e infantil na Índia

Fonte: PANDA, Basant Kumar; KUMAR, Gulshan; AWASTHI, Ashish. District level inequality in reproductive, maternal, neonatal and child health coverage in India. BMC public health, v. 20, n. 1, p. 1-10, 2020.

Modelos de Regressão Espacial

  • Hipótese de independência das observações em geral é Falsa \(\rightarrow\) Dependência Espacial

  • Efeitos Espaciais \(\rightarrow\) Se existir forte tendência ou correlação espacial, os resultados serão influenciados, apresentando associação estatística onde não existe (e vice-versa);

  • Como verificar ? \(\rightarrow\) Medir a autocorrelação espacial dos resíduos da regressão (Índice de Moran dos resíduos).

  • Autocorrelação espacial constatada ! E agora ? \(\rightarrow\) Utilizar Modelos de regressão que incorporam efeitos espaciais:

    1. Modelos Globais: utilizam um único parâmetro para capturar a estrutura de correlação espacial \(\rightarrow\) ex: Spatial Error Models (CAR)

\[y_i = \beta_0 + \sum^{p}_{k} \beta_k x_{ik} + \varepsilon_i \text{ , sendo } \varepsilon_i = \lambda W + \xi\]

  • Os efeitos da autocorrelação espacial são associados ao termo de erro. Sendo \(W\) a matriz de vizinhaça, \(\varepsilon\) a componente do erro com efeitos espaciais, \(λ\) é o coeficiente autoregressivo e \(\xi\) é a componente do erro com variância constante e não correlacionada.

  • A hipótese nula para a não existência de autocorrelação é que \(\lambda = 0\), ou seja, o termo de erro não é espacialmente correlacionado.

    1. Modelos Locais: parâmetros variam continuamente no espaço ex: Geographically Weighted Regression (GWR)

\[y_i = \beta_{0}(u_i, v_i) + \sum^{p}_{k} \beta_k (u_i, v_i) x_{ik} + \varepsilon_i \text{ , sendo } (u_i, v_i) \text{ as coordenadas geográficas}\] Fonte: ARDILLY, Pascal et al. Manuel d’analyse spatiale. 2018.

Serão feitas duas aplicações de Estatística Espacial utilização o pacote estatístico R:

Aplicação I: Utilizando a biblioteca tmap para construção de mapas temáticos

library(tmap)
data(World)
tmap_style("classic")
# Desenhando um rápido mapa temático mundial para esperança de vida.
qtm(World, fill = "life_exp")

Aplicação II: Baixando e construindo mapas a partir da biblioteca geobr

  • Bibliotecas que serão utilizadas:
library(ggplot2)
library(dplyr)
library(viridis)
library(geobr)
library(sf)
library(maptools)
library(leaflet)
library(ggspatial)
  • Para acessar os dados dos limites territoriais de todos os estados brasileiros é necessário utilizar a função read_state.
brasil.ufs <- read_state(code_state = "all", year=2019, showProgress = FALSE)
  • Primeiro, vamos fazer um gráfico apenas com as geometrias.
ggplot(brasil.ufs) + 
  geom_sf()

  • Para a construção de um mapa onde cada estado recebe uma cor de acordo com a sua região geogrática, procedemos da seguinte forma:
ggplot(brasil.ufs) + 
  geom_sf(aes(fill = name_region)) + 
  labs(fill="Região")

  • Agora utilizaremos dados relativos a porcentagem de municípios com rede de esgoto de acordo com a unidade da federação (fonte dos dados ) para fazer um gráfico. Vamos associar esses dados a tabela de acordo com a variável State e padronizaremos a porcentagem para variar de 0 a 100.
# Entrando com os dados observados na wikipedia
acesso_san <- data.frame(code_state = c(12, 27, 16, 13, 29, 23, 53, 32, 52, 21, 51, 50, 31, 15, 
                                   25, 41, 26, 22, 33, 24, 43, 11, 14, 42, 35, 28, 17), 
                         com_rede = c(0.273, 0.412, 0.313, 0.177, 0.513, 0.696, 1.000, 0.974, 0.280, 0.065, 
                                      0.191, 0.449, 0.916, 0.063, 0.731, 0.421, 0.881, 0.045, 0.924, 0.353, 
                                      0.405, 0.096, 0.400, 0.352, 0.998, 0.347, 0.129))
# Construindo o mapa com ggplot
brasil.ufs |> 
  left_join(acesso_san, by = "code_state") |> 
  mutate(com_rede = 100*com_rede) |> 
  ggplot(aes(fill = com_rede), color = "black") +
    geom_sf() + 
    scale_fill_viridis(name = "Municípios com rede de esgoto (%)", direction = -1) + 
  xlab("Longitude") + 
  ylab("Latitude") +
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "br") +
    theme_minimal()

  • Uma forma alteranativa de apresentar esses mesmos dados se dá pela apresentação de círculos com raios proporcionais a porcentgem de municípios com rede de esgoto no centroide de cada geometria (nesse caso, UF).
# Juntando o banco com os dados + malha
coord_pontos <- brasil.ufs |> 
                  left_join(acesso_san, by = "code_state") |> 
                  mutate(com_rede = 100*com_rede) |> 
                  st_centroid()
# Construindo o mapa com ggplot
ggplot(brasil.ufs)+ 
  geom_sf() + 
  geom_sf(data = coord_pontos, aes(size = com_rede), col = "blue", alpha = .65,
          show.legend = "point") + 
  scale_size_continuous(name = "Municípios com rede de esgoto (%)") + 
  xlab("Longitude") + 
  ylab("Latitude") +
  annotation_scale(location = "bl") +
  annotation_north_arrow(location = "br") +
    theme_minimal()

  • Uma alternativa interativa para trabalhar com mapas é com a utilização do pacote leaflet
data.frame(st_coordinates(coord_pontos), 
           com_rede = coord_pontos$com_rede, 
           UF = coord_pontos$name_state) |>
  leaflet() |> 
    addTiles() |>
    addCircleMarkers(~ X, ~ Y,
                     label = ~ as.character(paste0(UF, ": ", com_rede, "%")),
                     labelOptions = labelOptions(textsize = "13px"),
                     radius = ~ com_rede/10,
                     fillOpacity = 0.5)

Créditos ao Thiago Mendonça

Aplicação II: Análise exploratória utilizando os dados de dengue em Dourados/MS

OBS: Os dados e as malhas geográficas utilizadas nessa apresentação, estão disponíveis no seguinte endereço: (link)

  • Biliotecas do R que serão utilizadas
library(readr)
library(tidyverse)
library(sf)
library(maptools)
library(spatstat)
library(tmap)
  • Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS
unzip('dados/pop2010.zip', exdir='dados')
pop2010 <- read_csv('dados/pop2010.csv')

unzip('malhas/Setor_UTM_SIRGAS.zip', exdir='malhas')
setor.sf <- read_sf('malhas/Setor_UTM_SIRGAS.shp', crs = 31981)

unzip('malhas/contorno.zip', exdir='malhas')
contorno.sf <- read_sf('malhas/contorno.shp', crs = 31981)

OBS: Um aspecto muito importante dos dados espaciais é o sistema de referência de coordenadas (CRS) que é usado. Por exemplo, uma localização de (140, 12) não é significativa se você sabe onde está a origem e se a coordenada x está a 140 metros, quilômetros ou talvez graus de distância dela (na direção x). (link)

  • Fazendo um join com as tabelas com os setores censitários + população
setor.sf <- setor.sf %>% mutate (idsetor = as.numeric(CD_GEOCODI)) %>% inner_join(pop2010,by='idsetor')
  • Lendo e plotando os casos de dengue georreferenciados em Dourados/MS
unzip('dados/dengue_dourados.zip', exdir='dados')
casos <- read_csv('dados/dengue_dourados.csv')
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)
ggplot(setor.sf) + 
  geom_sf(fill = 'white', color='black') +
  geom_sf(data=casos.sf, color='red',size=1) +
  theme_void()

  • Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS
unzip('dados/lixo_dourados.zip', exdir='dados')
lixo <- read_csv2('dados/lixo_dourados.csv')
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)
ggplot(setor.sf) + 
  geom_sf(fill = 'white', color = 'black') +
  geom_sf(data=lixo.sf,color='blue',size=1) +
  theme_void()

Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS

  • Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.
lixo2.sf <- st_intersection(contorno.sf, lixo.sf)
ggplot(setor.sf) + 
  geom_sf(fill = 'white', color = 'black') +
  geom_sf(data=lixo2.sf,color='blue',size=1) +
  theme_void()

  • Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)
ggplot(setor.sf) + 
  geom_sf(aes(fill=pop)) + 
  geom_sf(data=casos.sf,color='red',size=0.7, aes(colour = "Caso"), 
          show.legend = "point") +
  geom_sf(data=lixo2.sf,color='salmon',size=1, aes(colour = "Lixo"), 
          show.legend = "point") +
  scale_fill_distiller(palette ="PuBu", direction = 1) +
  scale_colour_manual(values = c("Caso" = "red", "Lixo" = "slamon")) +
  theme_minimal()

  • Agora serão construidos buffers de 500m de distância ao redor de cada ponto de coleta de lixo. Isso é interessante para verificar se os casos de dengue ocorrem em um raio de até 500m de distância dos pontos de coleta de lixo. Ou seja, a pergunta é, será que a distância dos pontos de coeta de lixo influenciam a ocorrência do caso de dengue ?

Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.

lixo_buffer <- st_buffer(lixo2.sf, 500) 

ggplot(setor.sf) + 
  geom_sf(aes(fill=pop)) + 
  geom_sf(data=lixo_buffer,color='gray', fill = "transparent", size=0.4) +
  geom_sf(data=casos.sf, color='red',size=0.7) +
  scale_fill_distiller(palette ="PuBu", direction = 1) +
  scale_colour_manual(values = c("Caso" = "red", "Lixo" = "slamon")) +
  theme_minimal()

  • Represntando os casos e o lixo de forma interativa
tm_shape(setor.sf) +  tm_borders("black") +
  tm_shape(casos.sf) + tm_dots("red") +
  tm_shape(lixo.sf) + tm_dots("green") +
  tm_shape(lixo_buffer) + tm_borders("blue") +
  tmap_mode("view")
  • Convertendo o dado de pontos (padrão pontual) para dados de área
## conta casos por setor 
casos.sf$contador <- 1 
casos <- setor.sf |>  
  st_join(casos.sf) |> 
  filter(CLASSI_FIN == 1) %>%  ## seleciona somente os casos confirmados
  group_by(ID) |> 
  summarise(casos = sum(contador))

st_geometry(casos) <- NULL  ## remove atributos de geometria

## numero de depositos de Lixo por setor 
lixo.sf$contador <- 1 

lixo <- setor.sf |>  
  st_join(lixo.sf) |> 
  group_by(ID) |> 
  summarise(lixo = sum(contador))

st_geometry(lixo) <- NULL ## remove atributos de geometria

# Inserindo as contagens dos casos e de pontos de coleta de lixo no atributo com a geometria.
setor.sf <- setor.sf |> 
  left_join(lixo,by = 'ID') |> 
  left_join(casos,by = 'ID')
  • Plotando o mapa temático dos casos por setor censitário
plot(setor.sf['casos'])

  • Plotando o mapa temático dos pontos de coleta de lixo por setor censitário
plot(setor.sf['lixo'])

  • Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário
setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0 # Transformando os missings em zero

summary(setor.sf$tx)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.736   4.065   4.311  56.399
  • Plotando a distribuição da incidência em Dourados/MS
library(wesanderson)
pal <- wes_palette("Moonrise3", 20, type = "continuous")

ggplot(setor.sf) + 
  geom_sf(aes(fill = tx), color = 'black') +
  scale_fill_gradientn(colours = pal) +
  ggtitle("Taxa de incidência de Dengue") + 
  theme_void()

  • Kernel por atributos

Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.

Primeiramente é necessário dissolver os polígonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA

dourados.contorno <- st_union(setor.sf)
plot(dourados.contorno)

dourados.w <- as.owin(st_geometry(dourados.contorno))
  • Extraindo os centróides dos polígonos em Dourados/MS
centroides <- st_centroid(st_geometry(setor.sf))

# Transformando em os centróides em formato sp
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c('X','Y')

plot(centroides.sp)

  • Colocando os pontos no formato sp
centroides.ppp <- ppp(centroides.sp$X,centroides.sp$Y, dourados.w)

plot(centroides.ppp,pch=19,cex=0.5)

  • Fazendo o kernel por atributo da taxa de detecção
kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx, scalekernel = TRUE)
plot(kernel.tx)

  • Construindo a matriz de vizinhança para verificar a autocorrelação espacial
library(spdep)
viz <- poly2nb(setor.sf)
viz 
Neighbour list object:
Number of regions: 284 
Number of nonzero links: 1726 
Percentage nonzero weights: 2.139952 
Average number of links: 6.077465
  • Iremos precisar da coordenadas dos centróides
setor.sp <- as(setor.sf, 'Spatial') # convertendo em formato sp
coord <- coordinates(setor.sp) # coordenadas dos centroidas dos poligonos de dourados
class(setor.sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
  • Verificando a malha de conectividade da vizinhança de Dourados/MS
viz.sf <- as(nb2lines(viz, coords = coord), 'sf')
viz.sf <- st_set_crs(viz.sf, st_crs(setor.sf))

# Plota o grafo de conectividade por contiguidade
mapa.viz <- ggplot(setor.sf) + 
  geom_sf(fill = 'salmon', color = 'white') +
  geom_sf(data = viz.sf) +
  theme_minimal() +
  ggtitle("Vizinhança por \n conectividade") +
  ylab("Latitude") +
  xlab("Longitude")
mapa.viz

  • Obtendo a correlação da taxa de incidência de dengue Dourados/MS
pesos.viz <- nb2listw(viz)
moran.test(setor.sf$tx, pesos.viz)

    Moran I test under randomisation

data:  setor.sf$tx  
weights: pesos.viz    

Moran I statistic standard deviate = 15.724, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.524720303      -0.003533569       0.001128660
  • Plotando o correlograma
correl <- sp.correlogram(viz, setor.sf$tx, order = 8, method = "I")
correl
Spatial correlogram for setor.sf$tx 
method: Moran's I
           estimate expectation    variance standard deviate Pr(I) two sided
1 (284)  0.52472030 -0.00353357  0.00112866          15.7239       < 2.2e-16
2 (284)  0.23776816 -0.00353357  0.00048540          10.9525       < 2.2e-16
3 (284)  0.09536593 -0.00353357  0.00028358           5.8729       4.282e-09
4 (284) -0.02063378 -0.00353357  0.00019736          -1.2172         0.22351
5 (284) -0.07589158 -0.00353357  0.00016732          -5.5939       2.221e-08
6 (284) -0.06448448 -0.00353357  0.00016165          -4.7940       1.635e-06
7 (284) -0.02708340 -0.00353357  0.00016725          -1.8210         0.06861
8 (284) -0.02744543 -0.00353357  0.00018328          -1.7662         0.07735
           
1 (284) ***
2 (284) ***
3 (284) ***
4 (284)    
5 (284) ***
6 (284) ***
7 (284) .  
8 (284) .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(correl)

  • Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.
setor.sf$pval <- localmoran(setor.sf$tx, pesos.viz)[,5]

tm_shape(setor.sf) + 
  tm_polygons(col='pval', title="p-valores", breaks=c(0, 0.01, 0.05, 0.10, 1), border.col = "white", palette="-Oranges") +
  tm_scale_bar(width = 0.15)
  • Moran Local (Lisa Map) da taxa de incidência Dourados/MS
resI <- localmoran.sad(lm(setor.sf$tx ~ 1), 1:length(viz), viz, style = "W")
summary(resI)[1:10,]
    Local Morans I Stand. dev. (N)   Pr. (N) Saddlepoint Pr. (Sad)  Expectation
1 1   -0.009885292     -0.01575255 1.0125682 -0.03600714 0.9712767 -0.003533569
2 2    0.226981101      0.46511580 0.6418485  0.70056077 0.4835772 -0.003533569
3 3   -0.010910944     -0.01829621 1.0145974 -0.04041673 0.9677609 -0.003533569
4 4    0.484675519      1.31014141 0.1901480  1.43930439 0.1500643 -0.003533569
5 5    0.018648578      0.05952726 0.9525322  0.09384037 0.9252360 -0.003533569
     Variance    Skewness Kurtosis   Minimum  Maximum         omega       sad.r
1 1 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001369880 -0.01569409
2 2 0.2456263 -0.04188294 8.752890 -70.87400 69.87400  0.0028119325  0.44472643
3 3 0.1625854 -0.05147857 8.753481 -57.75455 56.75455 -0.0001590844 -0.01822753
4 4 0.1388594 -0.05570264 8.753780 -53.41199 52.41199  0.0066304485  1.08297547
5 5 0.1388594 -0.05570264 8.753780 -53.41199 52.41199  0.0005595065  0.05929966
          sad.u
1 1 -0.01569909
2 2  0.49831660
3 3 -0.01823490
4 4  1.59298211
5 5  0.05942125
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]
setor.sf$MoranLocal <- summary(resI)[,1] 

library(scales)

ggplot(setor.sf) + 
  geom_sf(aes(fill = MoranLocal), color = 'black') +
  scale_fill_gradientn(colours=c("blue", "white", "red"), 
                       values=rescale(c(min(setor.sf$MoranLocal), 0, max(setor.sf$MoranLocal))), guide="colorbar") + 
  ggtitle("Moran local") + 
  theme_void()

Bibliografia sugerida

  • Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds). Análise Espacial de Dados Geográficos. Brasília, EMBRAPA, 2004.

  • Interactive Spatial Data Analysis by Trevor C. Bailey , Anthony C. Gatrell Routledge, 1995

  • Applied Spatial Statistics for Public Health Data; Lance A. Waller, Carol A. Gotway Wiley-Interscience 1St ed. 2004

  • Applied Spatial Data Analysis with R; Roger S. Bivand, Edzer Pebesma , Virgilio Gomez-Rubio Springer; Edição: 2nd ed. 2013

  • Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow. Geocomputation with R, 2021.

Programa de Pós-Graduação em Epidemiologia em Saúde Pública (ENSP/FIOCRUZ)

  • Programa de Pós-Graduação em Epidemiologia em Saúde Pública (ENSP/FIOCRUZ) link)

Programa de Pós-Graduação em Ciências Veterinárias (IV/UFRRJ)

  • Programa de Pós-Graduação em Ciências Veterinárias (IV/UFRRJ) link

Obrigado !!!!!